#ifndef _MBPT_
#define _MBPT_
#include "../DF_Common/basis_state.h"
#include "../DF_Common/const_nucl.h"
#include "../DF_Common/coupling_coeff.h"
#include "../DF_MB/normal_order.h"
#include "../DF_Common/time_count.h"
//#include "operator.h"
#include <cmath>
#include <fstream>
#include <iostream>
#include <string>
#include <vector>
#include <numeric>
using std::string;
using std::vector;

class MBPT_SEC {
public:
  MBPT_SEC(SingleP &SP_x, DoubleP &DP_x, Norm_Ord_Control &norm_ct_x)
      : SP(SP_x), DP(DP_x), norm_ct(norm_ct_x) {}
  SingleP &SP;
  DoubleP &DP;
  DoubleP DP_Model; // Only Model space
  Norm_Ord_Control &norm_ct;
  double Energy_shift; // for U and U_HO defination of single particle state

  SixJ_coeff sixJ_x;
  vector<double> Sqrt_vec;

  double hbar_omega;
  int E_cut;

  int Derivative_Max = 10;
  double Norm_fac_C = std::sqrt(2);
  double EPS = 1e-8;
  int Sm_build_flag=-1;
  // one-body val & val_d derivative
  vector<double> vec_1B, vec_OA, vec_OB, vec_OC, vec_OD, vec_OE;
  vector<vector<double>> Sm_vec; // Sm = S/(m!)
  vector<double> Frac_vec;
  // two-body val & val_d derivative
  vector<vector<vector<vector<double>>>> Qm_vec, Sm2Bf1B_vec, vec_TwA_d, vec_TwB_d,
      vec_TwC_d, vec_TwD_d, vec_TwE_d, vec_TwF_d;

  int sp_core, sp_P_beg, sp_P_end, sp_Q_beg, sp_Q_end;
  int sp_core_Tz, sp_P_beg_Tz, sp_P_end_Tz, sp_Q_beg_Tz, sp_Q_end_Tz;

  void init();
  void build_First();
  void build_Second();
  // for check
  void print_U(string file);
  // one-body diagrams
  void build_sec_Sm();
  void print_sec_OB_sep();
  void print_sec_OB(string file);
  void print_sec_OB_deriva();
  void print_sec_OB_Luigi(string file);
  double cal_OneBody(int a, vector<double> &vec_d);
  // two-body diagrams
  void build_sec_Qm();
  void build_sec_TwB_sep();
  void print_sec_TwB();
  void print_sec_TwB_deriva();
  double cal_TwoBody(int a, int b, int c, int d, int J2, int Flag,
                     vector<double> &vec_d);
  double cal_TwoBody_normalized_Tz(int a, int b, int c, int d, int J2, int Flag,
                                   vector<double> &vec_d);
  double cal_2Bfr1B_normalized_Tz(int a, int b, int c, int d, int J2,
                         vector<double> &V_2f1B);
  // X-coefficient
  double cal_X0(int j2a, int j2b, int J2ab, int j2c, int j2d, int J2cd,
                int J2ac, int J2bd);
  // cross-coupled ME <ab|V|cd>
  // 1: Jc->Ja  Jd->Jb
  double cal_CRC_ME_I(int a_Tz, int b_Tz, int c_Tz, int d_Tz, int J2, int T2z);
  // 2: Jc->Jb  Jd->Ja
  double cal_CRC_ME_II(int a_Tz, int b_Tz, int c_Tz, int d_Tz, int J2, int T2z);
  double get2B_normalized_Tz(int a_Tz, int b_Tz, int c_Tz, int d_Tz, int J2,
                             int T2z) {
    double tbme = this->get2B(a_Tz, b_Tz, c_Tz, d_Tz, J2, T2z);
    double Norm_fac = 1.0;
    double Norm_fac_t = 1.0;
    if (a_Tz == b_Tz)
      Norm_fac_t = Norm_fac_t * (std::sqrt(2));
    if (c_Tz == d_Tz)
      Norm_fac_t = Norm_fac_t * (std::sqrt(2));
    Norm_fac = 1.0 / (Norm_fac_t);
    double tbme_norm = tbme * Norm_fac;
    return tbme_norm;
  }

  double get1B(int a_Tz, int b_Tz) { return norm_ct.get1B(a_Tz, b_Tz); }
  double get_U_HO(int a_Tz, int b_Tz) { return norm_ct.get_U_HO(a_Tz, b_Tz); }
  double get2B(int a_Tz, int b_Tz, int c_Tz, int d_Tz, int J2, int T2z) {
    return norm_ct.get2B(a_Tz, b_Tz, c_Tz, d_Tz, J2, T2z);
  }

private:
  // get1B

  // one-body diagrams
  double cal_OA(int a, double E_start, vector<double> &vec_d);
  double cal_OB(int a, double E_start, vector<double> &vec_d);
  double cal_OC(int a, double E_start, vector<double> &vec_d);
  double cal_OD(int a, double E_start, vector<double> &vec_d);
  double cal_OE(int a, double E_start, vector<double> &vec_d);
  // two-body diagrams
  double cal_TwA(int a, int b, int c, int d, int J2, int T2z, int Par,
                 double E_start, vector<double> &vec_d);
  double cal_TwB_1(int a, int b, int c, int d, int J2, int T2z, int Par,
                   double E_start, vector<double> &vec_d);
  double cal_TwB_2(int a, int b, int c, int d, int J2, int T2z, int Par,
                   double E_start, vector<double> &vec_d);
  double cal_TwC_1(int a, int b, int c, int d, int J2, int T2z, int Par,
                   double E_start, vector<double> &vec_d);
  double cal_TwC_2(int a, int b, int c, int d, int J2, int T2z, int Par,
                   double E_start, vector<double> &vec_d);
  double cal_TwD(int a, int b, int c, int d, int J2, int T2z, int Par,
                 double E_start, vector<double> &vec_d);
  double cal_TwE(int a, int b, int c, int d, int J2, int T2z, int Par,
                 double E_start, vector<double> &vec_d);
  double cal_TwF(int a, int b, int c, int d, int J2, int T2z, int Par,
                 double E_start, vector<double> &vec_d);
  //
  int Tri(int &j2a, int &j2b, int &j2c) {
    // I+J>=K, I+K>=J, J+K>=I,
    // I/2+J/2+K/2 = INTEGER.
    // TRI=1, WHEN TRIADIC CONDITION IS FULFILLED, TRI=-1 OTHERWISE
    int L2 = j2a + j2b + j2c;
    if (int(L2) != (L2 / 2) * 2)
      return -1;
    L2 = L2 / 2;
    if (j2a * j2b * j2c < 0)
      return -1;
    if ((L2 - j2a) * (L2 - j2b) * (L2 - j2c) < 0)
      return -1;
    return 1;
  }
  int phase_cal(int &n) const { return n % 2 ? -1 : 1; }
};

class EFF_KK {
public:
  EFF_KK(MBPT_SEC &MB_sec_t, SingleP &SP_t, DoubleP &DP_Model_t)
      : MB_sec(MB_sec_t), SP(SP_t), DP_Model(DP_Model_t) {}
  MBPT_SEC &MB_sec;
  DoubleP &DP_Model; // Only Model space
  SingleP &SP;
  double EPS = 1e-8;
  vector<double> Frac_vec;

  int Derivative_Max;
  int sp_core, sp_P_beg, sp_P_end, sp_Q_beg, sp_Q_end;
  int sp_core_Tz, sp_P_beg_Tz, sp_P_end_Tz, sp_Q_beg_Tz, sp_Q_end_Tz;

  vector<vector<vector<double>>> vec_2B_eff;
  vector<double> vec_1B_eff;

  vector<vector<vector<int>>> info_fold_vec;// F_i = { x terms }: [i][x][m] = the number of Q_m

  void init();
  void print();
  void build_info_fold();
  double H_next_cal(vector<double> &Qm_vec, double H_last);
  double H_iter(vector<double> &Qm_vec);
  double H_Fold(vector<double> &Qm_vec);
  double Heff_cal_2B(int ch, int index_F, int index_I);
  double Heff_cal_1B(int a);
  void build_2B();
  void build_1B();

private:
  int Tri(int &j2a, int &j2b, int &j2c) {
    // I+J>=K, I+K>=J, J+K>=I,
    // I/2+J/2+K/2 = INTEGER.
    // TRI=1, WHEN TRIADIC CONDITION IS FULFILLED, TRI=-1 OTHERWISE
    int L2 = j2a + j2b + j2c;
    if (int(L2) != (L2 / 2) * 2)
      return -1;
    L2 = L2 / 2;
    if (j2a * j2b * j2c < 0)
      return -1;
    if ((L2 - j2a) * (L2 - j2b) * (L2 - j2c) < 0)
      return -1;
    return 1;
  }
};

#endif
